LANDEN TRANSFORMATIONS AND THE INTEGRATION OF 
RATIONAL FUNCTIONS 
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Abstract. We present a rational version of the classical Landen transfor- 
mation for elliptic integrals. This is employed to obtain explicit closed-form 
expressions for a large class of integrals of even rational functions and to de- 
velop an algorithm for numerical integration of these functions. 



1. Introduction 
We consider the space of even rational functions of degree 2p 



l 2p := \R{z) = ^-\p{z)-.= Y. h ^ {P ~ 1 ' k) andQ(*):=E 
I ^^ z > fc=0 fe=0 



with positive real coefficients a/. , bu G K+ normalized by the condition do = a p = 1 , 
the space 

oo 

<£oo := |J e 2p 
P =i 

of normalized even rational functions, and the set of 2p — f parameters 
^2p := {a>i, • • • , flp-i; bo,- - ■ , bp-i}. 

We describe an algorithm to determine, as a function of the parameter set ^2p, a 
closed-form expression of the integral 



(f.l) I := / R{z)dz 

Jo 

for a large class of functions R £ £oo- The function R is called symmetric if its 
denominator Q satisfies Q(l/z) = z~ 2p Q(z). This is equivalent to its coefficients 
being palindromic, i.e. a,j = a p _j for 1 < j < p. 

The class of symmetric functions plays a crucial role in this algorithm. Define 

<£ s 2p := {R € <E2p den(R) is symmetric} 

(where den(R) denotes the denominator of i?), the class of rational functions with 
symmetric denominators of degree 2p, and 



p=i 
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For ragff define 

&2 P ■= {R G <Soo (den(R)) 1/{m+1) is even, symmetric of degree 2p} 

-2p S := ^2p ^ ^2p; so a function i? 6 ^^,' s can be written in the form 

P(z) 



R(z) = 



Q m+1 { z y 



where P(z) is an even polynomial and Q(z) is an even symmetric polynomial of 
degree 2p. 

The method of partial fractions gives (in principle) the value of I in terms of 
the roots of Q. Symbolic computations yield either a closed-form answer, an ex- 
pression in terms of the roots of an associated polynomial, or the integral returned 
unevaluated. 

The algorithm described here allows only algebraic operations on elementary 
functions and changes of variables of the same type. In particular, we exclude the 
solution of algebraic equations of degree higher than 2. We say that a rational 
function R £ is computable if its integral can be evaluated by our algorithm. 

The first step in the algorithm is to consider symmetric rational functions. In 
Section [2] we prove a reduction formula, i.e. a map $ p : £| p — > D\ p that reduces the 
computability of the integral of the symmetric function R to that of one of degree 
Tjdeg(R). Here 9\ p is the space of rational functions with denominator of degree p. 
These new functions are not necessarily symmetric, and this is the main limitation 
of our algorithm. The details of 3" p require the evaluation of some binomial sums 
which are presented in Appendix [X] The classical Wallis' formula 

dz 7r / 2m 

(z 2 + l) m+1 ~ 2 2m+1 V m 

shows that every R £ (£™ is computable. The reduction formula now implies that 
every R £ €™ is computable. This is described in Section [3] The computability of 
R G €4 is also a consequence of the classical theory of hypergeometric functions; 
the details are given in [3j . In Section 2] we compute the integral of every function 
in (£g l,s , where the reduction method expresses these integrals in terms of functions 
in (£™. The algorithm does not, in general, provide a value for the integral of a 
nonsymmetric function of degree 8. 

The final piece of the algorithm is described in Section [5] for R € (£2 P , the 
symmetrization of its denominator produces a (symmetric) rational function in (£.4 P 
with the same integral as R. The reduction formula now yields a new function in 
<&2p with the same integral as R. We thus obtain a map T2 P : £ 2 p — * ^2 P such that 

/•oo />oo 

(1.2) / R{z)dz = I 1 2p {R{z))dz. 

Jo Jo 

In particular, the class of computable rational functions of degree 2p is invariant 
under forward and backward iteration of S^p- This map is the rational analog of 
the original Landen transformation for elliptic integrals described in [H [13j . The 
map 1 2p can also be interpreted as a map on the coefficients $ 2 p : &2p — * ^2p 
where D£ p = x R + p . 
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The case of $6 is described in detail in Section [6] and is given explicitly by 

9 + 5ai + 5a 2 + a\a 2 



(1.3) oi 

«2 
&0 

h 

l>2 



(oi + a 2 + 2) 4 /3 
oi + a2 + 6 
(ai + a 2 + 2) 2 / 3 

&0 + &1 + ^2 

(ai + a 2 + 2) 2 / 3 
6o(o 2 + 2) + 2^+62(01 +3) 
01 + a 2 + 2 
&o + h 



ii+a 2 + 2) 1 /3- 



Let xo := (oi, 02; bo, &i, 62)- Then $6 : 0@ — > Dg is iterated to produce a sequence 
x n+ i := $ 6 (x„) of points in that yield a sequence of rational functions with 
constant integral. We have proved in [3] the existence of L £ R + , depending upon 
the initial point Xo = (ax, a 2 ; bo, b\, b 2 ), such that x„ — * (3, 3; L, 2L, L). Thus 

n f°° b Q z i + b l z 2 + b 2 ir 
(1.4 / —s -. s dz = L xn x — . 

V ' J z e + ai z 4 + a 2 z 2 + 1 V " J 2 

This establishes a numerical method to compute the integral in (jl.4|) . 

Numerical studies on integrals of even degree 2p suggest the existence of a lim- 
iting value L = £( x o) such that the sequence x„ := $2p(x, l ) satisfies 

p\ fp\ ( v \ j P -i\ ( P -i\ ( P -i 



V V 2 / vp - v \ / V 1 / V» - 1 

The integral of the original rational function is thus J x £(xo). The proof of 
convergence remains open for p > 4. Examples are given in Section [7J 

The most important issues left unresolved in this paper are the convergence of 
the iteration of the map <& 2p discussed above and the geometric interpretation of the 
Landen transformation 1 2p . Finally, the question of integration of odd functions 
has not been addressed at all. 



Some history. The problem of integration of rational functions R(z) = P(z) /Q(z) 
was considered by J. Bernoulli in the 18 t/l century. He completed the original 
attempt by Leibniz of a general partial fraction decomposition of R(z). The main 
difficulty associated with this procedure is to obtain a complete factorization of 
Q(z) over K. Once this is known the partial fraction decomposition of R(z) can be 
computed. The fact is that the primitive of a rational function is always elementary: 
it consists of a new rational function (its rational part) and the logarithm of a second 
rational function (its transcendental part) . In his classic monograph [9 G. H. Hardy 
states: The solution of the problem (of definite integration) in the case of rational 
functions may therefore be said to be complete; for the difficulty with regard to the 
explicit solution of algebraical equations is one not of inadequate knowledge but of 
proved impossibility. He goes on to add: It appears from the preceding paragraphs 
that we can always find the rational part of the integral, and can find the complete 
integral if we can find the roots of Q{z) = 0. 



4 



GEORGE BOROS AND VICTOR H. MOLL 



In the middle of the last century Hermite [10] and Ostrogradsky [15] developed 
algorithms to compute the rational part of the primitive of R(z) without factor- 
ing Q(z). More recently Horowitz [TT] rediscovered this method and discussed its 
complexity. The problem of computing the transcendental part of the primitive 
was finally solved by Lazard and Rioboo [12] . Rothstein [17] and Trager [18]. For 
detailed descriptions and proofs of these algorithms the reader is referred to [5] and 
i- 

2. The reduction formula 

In this section we present a map $ p : (Sjp 5 — * &p that is the basis of the 
integration algorithm described in Section [5] The proof is elementary and the 
binomial sums discussed in the Appendix are employed. 

Let D p (z) be the general symmetric polynomial of degree 4p. We express the 
integral of z 2rl /£>™ +1 as a linear combination of integrals of z 2j /£'™ +1 where E p is 
a polynomial of degree 2p whose coefficients are determined by those of D p . 

Theorem 2.1. Let m,n,pe N. Define 

p 

(2.1) D p (d 1 ,d 2 ,--- ,d p ;z) = Y,d p+1 _ k (z 2k + z ip - 2k ) 

k=0 

and 

E p (di,d 2 ,--- ,d p ;z) = 




£ 2 2i-x z 2 iv -i) 3+i-l ( 3 + 2< - 2 j d 



»=i j=i 
(2.2) 

for d l S R+. Then for < n < (to + l)p - 1, 

z 2 " dz 



o ,d p ;z)) m+1 



U V 2j (E p (d u ---,d p ;z)) m+1 

(2.3) 

and for (to + l)p — 1 < n < 2p(m + 1) — 1 we employ the symmetry rule 

(2-4) N np = N 2p (jn+l)-l-n,p- 

Proof. First observe that (|2.4|) follows from the change of variable z — > 1/z. Now 
consider 

z 2 ™ dz 



(2.5) Nnvidi, - • • ,dr,:m) .— , 

/o (ELo^ +1 - fc (^- 2fc + ^)) m+ 
for < n < (m + l)p — 1. The substitution z = tan 8 yields 

»tt/2 _ £r2^n£r4(m+l)p-2n-2 ^ 







ELo d p+i-k {(1 - C 2 ) 2 f- fc C7 2fe + (1 - C 2 ) fe C7 4 P- 2fe } ) r 
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where C — cos 8. Letting tp — 29 and D = cost/' = 2C 2 — 1 then gives 

(1 - D 2 ) n {\ + £>)2(m+l)p-2n-l ^ 



N 

1 y n,p 



'° ( ELo dp+i-fc(l - D 2 ) k {(1 - Dfp-^ + (1 + Dfp-^} 
Now observe that the integrals of the odd powers of cosine vanish when we expand 

(1 + £))2(m+l)p-2n-l ) producing 



N nn = 2 



— in 



tt/2 (1 _ £)2-)n ^(m+l)p-7i-l (2( m +l)p-2n-l^ D 2 3 



{n=o^fi-*(i - ^ 2 ) fc E^o cy*)s*} m+1 ' 

A second double angle substitution ip = 2ip gives 

(1 - £)" ^2(m+l)p-2n-l^ 2 p(m+l)-n-i-l^ + 



{ELo <W-*(1 - E) k (V ) 2P_k_i ( 1 + M 

where £7 = costp = 2D 2 — 1. The change of variable z — tan(<y9/2) then yields 

„oo z 2n yv(m+l)p-n-l /2(m+l)p-2n-l\ _|_ ^2\(m+l)p-n-j-l ^ 
_ 2 — m / ^3=0 v _ . 

Jo {ELo^ + i-^ 2fe (E^o(V' s )( 1 + ^) p ^)} W+1 

(2.6) 

Finally, we modify (|2.6p using Lemma [A.2l and Lemma lA.4l with N — (m+l)p—n—l 
to produce EM- □ 



Note that the previous theorem associates to each rational function of symmetric 
denominator 

R(z] b s z 2s + bs-iz**- 1 ) + ■ ■ ■ + b 

(z i P + d p z i P- 2 + --- + 2d 1 z 2 P + --- + l) m+1 

a new rational function 

(m+l)p-l (m+l)p-n-l o-n-1 + A ^2((m+l)p-l-i) 

h U ^ 2j J(E p (d 1 ,---,d pi z)) m+1 

such that 



/ R!(z)dz = / R 2 (z)dz 
Jo Jo 



3. The quartic case 

In this section we describe the computability of rational functions R G €™. 
These are functions of the form 

P(z) 

^ _ (z 4 + 2az 2 + 1)™+! 

where P(z) is an even polynomial of degree 4to + 2. Observe that the normalization 
o-a = 0-2 = 1 makes the denominator of R automatically symmetric. It suffices to 
evaluate 

(3.1) N nA (di;m) := 



(z 4 + 2d 1 z 2 + l)"»+i 
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where < n < 2m+l is required for convergence. From (12. 4|) we have N n ^(di; m) = 
N2m-i-n,i{d\; m), so we may assume < n < m. We now employ Theorem 12. II to 
obtain a closed form expression for N n ^{d\;m). 

Theorem 3.1. Let m G N and assume < n < m. Then 

f°° z 2n dz 

(3.2) N^dum) 



o (z 4 + 2d lZ 2 + 1)™+! 

2 3m+3/2( 1 + a )m+l/2 X 2 J (1 



j=0 

For m + 1 < n < 2m + 1 we have 

z 2n dz 

(3.3) ' 



2m - 2J - 1 \ /m - n + j \ / 2j\ / m 



o (z 4 + 2d!Z 2 + l) m+1 



7T ^ w \j ^ 2m — 2j — 1\ /m — n + A ( 2 A / m\ 

23».+3/ 2( l + dl)m+ l/2 X A. 2 ^ 1 + ^ X (, m -j ){ 2 j 

Proof. We apply the result of the Theorem 12 . II with D\(di \ z) = z 4 + 2d\z 2 + 1 and 
Ex{d\\z) = (1 + dx)z 2 + 2, so that 

f°° z 2 " dz f m - n +j\ f°° z 2{ - m -^dz 

J (z 4 +2d l2 2 + l)m+l [ 2 j )J 



3=0 

The change of variable u = (1 + d\)z 2 /2 then yields 



o + + 2) 



m+l ' 



y 2(m— j) j r /■oc m —j—i/2 j 

2 -U+3/2)^ 1 + di y m +j-l/2 I " ''" 



o ( (l + di)z 2 + 2 ) m + 4 " v y (l + u)™+! 



^2m - 2 A /2j^ /m^ 1 2 _(2„ l+3 +3/2) ( ^ 1 + di y(m-j+l/2) ^ 



where we have used 



00 if- 1 ' 2 , 7T /2r\ (2{s - r - 1)\ (s - 1\ 1 

au — iii ii 



q (l + u) s 2 2 ( S - 1 > V r / V 8-r-l 



The algorithm also requires a scaled version of iVo^di; m). 

Corollary 3.2. Let b > 0, c > 0, a > -*\/&e, m € N, and < n < m. Define 

z 2n dz 



□ 



N nA (a, b, c;m) 
Then for < n < m, 



(fez 4 + 2az 2 + c) 
-1/2 



m+l ' 



N nA (a,b,c;m) = tt ^c(c/6) m -" |8(a + v 7 ^)} 2 ™^ 



1. /2m - 2fc\ fm~n + k\ ( 2h\ I ' ra\ 1 / // 



X 

fc=0 

(3.4) 



\ rn — k J \ 2k J \ k J \ k J \ y/bc 
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and for m + 1 < n < 2m + 1, 

N nA (a,b,c;m) = n (c{c/b) m - n {8(a + Vbc")} 2 '"^ 



n x^\k /2m - 2k\ fm - n + k\ (2k\ fm\ ~ 1 f a x 



% m — k J \ 2k J \ k J \ k J \ *.fhf- 
fc=o v /v / \ / \ / w«i- 

(3.5) 

Proof. Let < n < m. The substitution u = z(6/c) 1 / 4 yields 
(3.6) 7V„, 4 (a,6,c;TO) = ——^^^— N^ 



so (|3.4[) then follows from Theorem 13.11 From (|2.4ll we have N ny 4(a,b,c;m) = 
N2m+i-n,4(c, b, a; m) for m + 1 < n < 2m + 1, giving (|3.5|) . □ 

4. The symmetric case of degree 8 

In this section we prove the computability of the set (£™' s of symmetric rational 
functions with denominator of degree 8 and establish an explicit formula for the 
integral 

f°° z 2n dz 

N nt s(ai,a 2 ;m) 



lo ( z& + a 2Z e + 2a x z A + a 2 z 2 + l) m+1 
where < n < Am + 3 is required for convergence. Observe that (|2.4p reduces the 
discussion to the case < n < 2m + 1. The expression (|2.2p . with p = 2, produces 
E 2 (ai,a 2 ; z) = (1 + a x + a 2 )z 4 + 2(a 2 + 4)z 2 + 8. 

Theorem 4.1. Every function in (£™' s is computable. More specifically, define 
ci := a 2 + 4, c 2 := 1 + a x + a 2 , and 

i W (m,nio 1)fls ) := 7r2-^+ 2+k +^ 2 c^- k ~ n/2 ( Cl + V^?-^ 1 ' 2 x 

/ 4to — n — k + 2\ ( 2m — 2j\ fm — k + j\ ( 2j\ / m\ 1 



k — n J \ m — j ) \ 2j J \ j J \j 
Then for m + 1 < n < 2m + 1. 1 + a\ + a 2 > and a 2 + 4 > -8^/8(1 + a x + a 2 ), 

„oo z 2n dz 2m+U-m-l 

/ (z8 + a2 z6 + 2aiZ 4 + a2Z 2 + 1)m+ i = E E twKnjai.aa), 
and for < n < m, 



(z 8 + a 2 z 6 + 2a 1 z 4 + a 2 z 2 + l) m +! 

m m — k 2ra+l k— m — 1 

^^{y(m,n;a 1) a 2 )+ (m, n; ai, a 2 ). 

fc=n J— fc— m+1 j— 

Proof. The reduction formula yields 



/o (z 8 + a 2 z 6 + 2aiz 4 + a 2 z 2 + l) m + 1 
(41) 23m+2 2 g 1 2 - 2 ^4m-n-fc + 2\ z 2fe 



lo (c 2 z 4 + 2c x z 2 + I 



m+1 
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We then use Corollary [32] to evaluate (l4~Tj) . □ 

5. A SEQUENCE OF LANDEN TRANSFORMATIONS 

The transformation theory of elliptic integrals was initiated by Landen in 1771. 
He proved the invariance of the function 

fV2 rlft 



G(a,b) , 

V» cos 2 8 + b 2 sin 6 
under the transformation 

(5.2) ai = (a + b)/2 h = Vab, 
i.e. that 

(5.3) G(ai,6i) = G(a,b). 

Gauss 7} rediscovered this invariance while numerically calculating the length of a 
lemniscate. An elegant proof of (|5.3[) is given by Newman in |14j . Here, the sub- 

stitution x = otan# converts 2G(a,b) into the integral of [(a 2 + x 2 )(b 2 + x 2 )] 
over R; the change of variable t = (x — ab/x)/2 then completes the proof. 

The Gauss-Landen transformation can be iterated to produce a double sequence 
(a n ,b n ) such that < a n — b n < 2~ n . It follows that a n and b n converge to 
a common limit, the so-called arithmetic- geometric mean of a and b, denoted by 
AGM(a, b). Passing to the limit in G(a, b) = G(a n , b n ) produces 

frr/2 

(5.4) 



2AGM(a,b) J ^/ a 2 cos 2 6 

+ b 2 sin 9 

The reader is referred to [1] and [13] for details. 

The goal of this section is to produce a map %2p ■ ^2p - ► &2p that preserves the 
integral, i.e. 



(5.5) / R{z)dz = / % 2 p{R{z))dz. 
Jo Jo 

This map is the rational analog of the original Landen transformation (|5.2[) . 

Theorem 5.1. Let R(z) = P(z)/Q(z) with 

p-i p 

(5.6) P(z) = J2 b i z2{p ^ 3) and Q( z ) = J2 a i z2(P ~ 3) - 

3=0 ]=0 

Define aj = for j > p, bj = for j > p — 1, 

j 

(5.7) dp+i-j = ^ a p -kUj-k 

k=0 



for < k <p- 1, 
(5.8) di = 



p 

2 



2 

A:=0 



a p—ki 



2p-l 

(5.9) Cj = a J b P-i-j+k 



fe=0 
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for < j < 2p — 1 , and also 

(5.10) Op(i) = 

Let 
(5.11) 
for 1 < i < p — 1 , and 

(5.12) b+ = Q{1) 2i 'p +1 'p- 2 



i + ELi 4 if * = o. 

Op(i) 



+ _ 



2 2i Q(l) 2 ( 1 - l /p) 



-l-i 
L fc=0 

for < i < p — 1. Finally, define the polynomials 



p — 1 — k + i 
2i 



p-i 



(5.13) 



p + w=E 6 



and Q + (z) = ^ a 



+ z 2(p- 4 )^ 



fc=0 



Then % 2p {R{z)) ■= P + (z)/Q+(z) satisfies fOjl . i.e. 



(5.14) 



P(z) 

o Q(*) 



50 P+(z) 



tfe. 



Proof. The first step is to convert the polynomial Q(z) to its symmetric form: 



I := 



P(z) 



o Q(z) 7 £>(*) 



C7(*) 



dz 



with 



2p-l 



C(z) - P{z) x z 2 *>Q(l/z) := J] 



Cfe -z 



_'/.■ 



fc=0 

p 



D(z) = Q(z) x z 2p Q(l/z) —^dp+i-k^- - - 



2fe , r 2(2p-fe) 



Then 



k=0 



roc z 2k dz 



E 

fc=0 



o Q s (z) 



Now employ the reduction formula in Section [2] to evaluate 

f°° z 2k dz 

Li 



lo Q s iz) 

Observe that one needs to evaluate Lk only for < k < p — 1. Indeed, the usual 
symmetry rule yields = Li v -\-k- The reduction formula now gives 



p-i-fc 

E ^ 

3=0 

p—k 



P — 1 — k + j 

2j 



z 2(p-l-j) dz 



i 



E 22( " 



i) 



p - + j - 2 
2j-2 



"p(p) 

with a p (i) as in (|5.10p and A = [a p (p)/a p (0)] 1 ^ 2p 



00 z 2 <P-J) dz 



□ 
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Note. The extension of this transformation to the case of 

P(z) 



)m+l( z ) 



dz 



requires explicit formulae for the coefficients of P(z) x ( z 2p Q(l/z)) m+1 and Q' m+1 (z) 



x 



?2p 



Q(V*)) 



m+l 



An algorithm for integration. Let x = (a, b) with a = (a\,--- ,a p _i), b 

>2p 



(bo, ■ ■ • , fep-i), and let D 2 = R+ p 1 x We then have a map 



x:=(a,b) -> x+:=(a + ,b+) 

where af and 6+ are given in (|5.111 I5.12[) . Iteration of this map, starting at x , 
produces a sequence x n+1 := $> 2p (x n ) of points in D 2p . The rational functions 
formed with these parameters have integrals that remain constant along this orbit. 
Numerical studies suggest the existence of a number L — L(x ) G such that 

Thus the integral of the original rational function is J x L. 



6. The sixth degree case 

We discuss the map 1 2p '■ £ 2 p — > <£ 2 p for the case p = 3. The effect of Tg on 
the coefficients ^p 6 = {^cb bi,b 2 , ai, a 2 } is denoted by $ 6 : 0^ — + Dg and is given 
explicitly by 



(6.1) at 

h 



9 + 5di + 5a 2 + aia 2 

(01 + a 2 + 2) 4 /3 

ai + a 2 + 6 
(ai + a 2 + 2) 2 /3 

&0 + &1 + &2 

(ai + a 2 + 2) 2 / 3 
fr (a 2 + 2) + 2fr 1 +6 2 (a 1 +3) 
ai + a 2 + 2 
&o + h 



2 



( ai +a 2 + 2)V3 
using Theorem 15. 11 The map $6 preserves the integral 

(6.2) U 6 (ai,a 2 , b ; b 1 ,b 2 ) ■= / -g- ■ 7- 2 —; dz 

J z° + a\z^ + a 2 x z + 1 

and the convergence of its iterations has been proved in [3] , the main result of which 
is the following theorem. 

Theorem 6.1. Let x := (a?, ag; &g, 6?, &g) S M+ 5 . Define x„+i := $ 6 (x„). Then 
C/6 is invariant under $6. Moreover, the sequence {(a", )} converges to (3, 3) and 
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{(&q , 6™, b 2 )} converges to (L, 2L, L), where the limit L is a function of the initial 
data xo. Therefore 

b z 4 + b 1 z 2 + b 2 . . tt 

dz = -L(xoJ x 



z b + aiz 4 + a 2 z 2 + 1 2 

This iteration is similar to Landcn's transformation for elliptic integrals that has 
been employed in [3] in the efficient calculation of n. Numerical data indicate that 
the convergence of x„ is quadratic. The proof of convergence is based on the fact 
that $ 6 cuts the distance from (ai, 02) to (3, 3) by at least half. 

A sequence of algebraic curves. The complete characterization of parameters 
(01,02) in the first quadrant that yield computable rational functions 



R(z) := 



b Q z 4 + 6iz 2 + b 2 



z 6 + a\z A + a 2 z 2 + 1 

of degree 6 remains open. The polynomial z 6 + a\z 4 + a 2 z 2 + 1 factors when 
ai = a 2 so the diagonal A := {(<2i, a 2 ) G M+ x R + : a± = a 2 } produces computable 
functions. In view of the invariance of the class of computable functions under 
iterations by $67 the curves X n := <P ( f n \A), with n G Z, are also computable. 
The curve Xi has equation 

(9 + 5ai + 5a 2 + aia 2 ) 3 = (ai + a 2 + 2) 2 (a 1 + a 2 + 6) 3 

and consists of two branches meeting at the cusp (3, 3). In terms of the coordinates 
x = a,i — 3 and y = a 2 — 3 the leading order term is Ti(x, y) — 1728(x — y) 2 . This 
curve is rational and can be parametrized by 

(6.3) ai{t) = t- 2 (t 5 -t 4 + 2t 3 -t 2 + t+l) 

a 2 (t) = t _3 (< 5 + t 4 - t 3 + 2t 2 - t + 1). 

The rationality of X n for n 5^ 1 and its significance for the integration algorithm 

remains open. The complexity of these curves increases with n. For example, the 

/ 2) 

curve X 2 := $3 (A) is of total degree 90 in x — a\ ~ 3 and y — a 2 — 3 with leading 
term 

T 2 (x,y) := 2 121 3 35 {x - y) 18 [-163(x 4 + y 4 ) + 668xy{x 2 + y 2 ) - \07Ax 2 y 2 } . 

The diagonal A can be replaced by a 2-parameter family of computable curves 
X(c, d) that are produced from the factorization of the sextic with a\ = c + d and 
a 2 = cd+ 1/d. All the images $g X(c, d) with neZ are computable curves. The 
question of whether these are all the computable parameters remains open. 

7. Examples 

In this section we present a variety of closed-form evaluations of integrals of ra- 
tional functions. 

Example 1. The integral 

23698523 tt 



dz — 



(z 4 + 4z 2 + l) 9 12230590464^6 
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is computed by Mathematica 3.0 using (|3.2|) in .01 seconds. The direct calculation 
took 12.27 seconds and 6.4 extra seconds to simplify the answer. 

Example 2. The integral of any even rational function with denominator a power 
of an even quartic polynomial can be computed directly by using Corollarv l3.2l For 
example: 

z 6 dz _ 11tt(14229567+ 4937288 y/6) 



(2z 4 + 2z 2 + 3) 11 440301256704 (1 + ) 21 / 2 

Example 3. The case n = in Q3.2| deserves special attention: 
( 7 - X ) N 0A (a;m) = 2m+3/2 ^J +1)m+1/2 P m (a) 



where 

,7.2, p. ( .) . ^g^frf)^:*)^!)' 

The polynomial P m {a) has been studied in pQ and [2]. 



Example 4. The case n = m in (|3.2|) yields 

_ f°° z 2m dz ir (2,n 

N mA {a;m) - (z 4 + 2az l + 1)m +l ~ 2 3m+3/2 (1 + a )m+l/2 x 

The change of variable z — > ^/z converts this integral to 

1 r°° z m-i/2 dz 

2 7 (z 2 + 2az+l)™+ 1 ' 

which is [S] 3.257.9. 



N mt i(a;m) = 



Example 5. A symmetric function of degree 6. The integral 



(x 6 + 4x 4 + 4x 2 + l) 5 ' J [(x 2 + l)(x 4 + 3x 2 + l)] 5 







can be computed by decomposing the integrand into partial fractions as 
1 1 6 11 31 



(x 2 + l) 5 (x 2 + l) 4 (x 2 + l) 3 (x 2 + l) 2 (x 2 + l) 
1 2x 2 4 3x 2 

+ 



(x 4 + 3x 2 + l) 5 (x 4 + 3x 2 + l) 5 (x 4 + 3x 2 + l) 4 (x 4 + 3x 2 + l) 4 
12 6x 2 32 14x 2 



(x 4 + 3x 2 + l) 3 (x 4 + 3x 2 + l) 3 (x 4 + 3x 2 + l) 2 (x 4 + 3x 2 + l) 2 

73 31x 2 



(x 4 + 3x 2 + l) (x 4 + 3x 2 + l)' 
Each of these terms is now computable yielding 

1407326 v / 5 - 3146875 

I = X 7T. 

160000 
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Example 6. Non-symmetric functions of degree 6. In this case we can use the 
scheme (|6.1j) to produce numerical approximations to the integral. For example, 
the evaluation of 

r °° 45z 4 + 25000z 2 



is shown below: 



1230 



z 6 + z 4 + 3000z 2 + 1 



dz 



n 


< 


a™ 


h n 
u 


V{ 







1 


3000 


45 


25000 


1230 


1 


.415786 


14.4465 


126.233 


63.2884 


88.3741 


2 


2.06562 


3.17262 


42.2607 


156.015 


83.6896 


3 


2.98142 


3.00338 


75.3541 


137.717 


65.1111 


4 


2.99999 


3. 


69.6338 


139.925 


70.2771 


5 


3. 


3. 


69.9589 


139.914 


69.9555 


6 


3. 


3. 


69.9572 


139.914 


69.9572 


7 


3. 


3. 


69.9572 


139.914 


69.9572 



Thus L ~ 69.9572 and 

45x 4 + 25000x 2 + 1230 



x 6 + x 4 + 3000x 2 + 1 



dx 



69.9572 x - ^ 109.889. 
2 



Example 7. Symmetric functions of degree 8. These integrals can be evaluated 
using Theorem 14. 11 For example: 

dz 



(14325195794 + 2815367209 V26 ) tt 



(z 8 + 5z 6 + Uz 4 + 5z 2 + I) 4 



14623232 (9 + 2V26) 7/2 



Example 8. As in the case of degree 6 we can provide numerical approximations 
to nonsymmetric integrals of degree 8. The iteration (|6.1[) is now replaced by 



n+l 



n+l 



un+1 
°0 



n+l 



un+1 
°2 



un+1 
°3 



4a™ a% 



10(a™ 



-8 (a" 



2) 



« + a" + a™ + 2) 3 / 2 
a> n + 6(a" + a") + 2(a' 2 l + 10) 



a" + a% + 2 



a" + a£ + I 



K 



J 2 " 

6™ 



2)V2 

-W 



2) 3 / 4 



b"(3a" + ag + 6) + b"(a" + 4) + b"(a" + 4) + 6"(3a" + a| + 6) 
K + + a£ + 2) 5 /4 

&3K + 5) + &' 2 l + + &ft(ag + 5) 



(a? 



2) 3 / 4 



« + a' 2 1 + 4 l + 2) 1 /4 
with initial conditions a°, a^i a 3! ^2' ^3- Then 

r b x e 

J x 8 + aix 6 



(7.3) U s (ai,a 2 ,a 3l b Q ,bi,b 2 ,b 3 ) 



hx 4 



b 2 x 2 



a 2 x^ 



a 3 x 2 + 1 



dx 
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is invariant under these transformations. 

Note. Numerical calculations show that (a™, a^, a^) — > (4, 6, 4) and that (6g , 6™, b%, b^) 
(1, 3, 3, 1) L for some L depending upon the initial conditions. 



Example 9. A symmetric function of degree 12. We use Theorem 12. II to evaluate 

z 18 dz 



I := 



(7.4) I = 



(z 12 + 14z 10 + 15z 8 + 4z 6 + 15z 4 + 14z 2 + l) s 
25tt(25V56- 54) 



301989888 



Here p = 3, n = 9, and m — 2, so n > (m + l)f> — 1 and we need to apply the 
transformation z —* 1/z to reduce the value of n. Indeed, we have 



I = 



z 16 (iz 



(z 12 + 14z 10 + 15z 8 + 4z 6 + 15z 4 + 14z 2 + l) 3 : 
and Theorem 12.11 now yields 

z 16 dz 



I 



-17 



131072(1 + z 2 ) 3 (l + 4z 2 + z 4 ) 3 ' 



The new integrand is expanded in partial fractions in the variable t = z 2 to produce 

El. 



Example 10. We use Theorem 12. II to evaluate 

/ := 

where 



o Q 2 (z) 



Q{z) = z 20 + 6z 18 + 93z 16 - 24z 14 + 162z 12 + 548z 10 + 162z 8 - 24z 6 + 93z 4 + 6z 2 + 1. 
The factorization 

Q(z) = (l + z 2 ) 2 T(z)T(-z) 

with 

T(z) = z 8 - 2z 7 + 4z 6 + 14z 5 + 6z 4 - 14z 3 + 4z 2 + 2z + 1 

leads to a partial fraction expansion containing the term 

72 - 501z + 1994z 2 - 2617z 3 + 1228z 4 - 43z 5 + 34z 6 - 55z 7 
8388608(1 - 2z + 4z 2 + 14z 3 + 6z 4 - 14z 5 + 4z 6 + 2z 7 + z 8 ) 2 ' 

which we were unable to integrate; furthemore, the roots of T(z) = cannot be 
evaluated by radicals. The procedure described in Theorem 12.1) however, shows 
that 



I = 



z 10 (4 + z 2 )(z 6 + 36z 4 + 96z 2 + 64) dz 
524288(z 2 + l) 2 (z» + 3z6 + 8z 4 + 3z 2 + l) 2 ' 
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the integrand of which can be expanded to yield 

9 f 00 dz 75 



I 



dz 



8388608 J (z 2 + l) 2 8388608 

1921z 6 + 10815z 4 + 4111z 2 + 1462 
o 2097152 (z 8 + 3z 6 + 8z 4 + 3z 2 + l) 2 

00 91z 6 + 719z 4 + 1259z 2 - 5764 



z 2 + l 



dz 



dz 

o oooouuo(z 8 + 3z 6 + 8z 4 + 3z 2 + l) 
Every piece is now computable, with the final result 

_ (6480 - 509V15) tt 
24159191040 

Example 11. The symmetric functions of degree 16 have denominator 
D i {d 1 ,d 2 ,d 3 ,d i ;z) = z 16 + d 4 z 1A + d 3 z 12 + d 2 z 10 + 2d x z* + d 2 z 6 + d 3 z 4 + d 4 
the integral of which is computed in terms of 

£■4(^1,^2,^3,^4; 2) = (1 + di + d 2 + d 3 + d 4 )z 8 + 2(16 + d 2 + 4d 3 + 9d 4 )z 6 
+ 8(20 + d 3 + 6di)z 4 + 32(8 + d 4 )z 2 + 128. 
This new integral is symmetric provided 



(7.5) 

Introduce the new parameters 



15 
112 



3 

-4 



d 4 . 



e j — ">j 
Then (f7T5|) yields 
(7.6) 



5-j 



for 2 < j < 4 and ei = d\ - \ 



ei 




' 3 " 




"-8" 


e 2 _ 




-4 


e 3 + 


7 



P4- 



Thus, if the original denominator has the form 

D A (z) = (z 16 + 1) + d 4 (z 14 + z 2 ) + d 3 (z 12 + z 4 ) + (112 - 4d 3 + 7 d A ){z w + z 6 ) 
+2(15 + 3c? 3 - 8d 4 )z s , 

the integral 

P(z) 



\ rn+1 



dz 



h (D 4 (z)Y 

is reduced to an integral with symmetric denominator of degree 8 and these are 
computable. We can thus evaluate a 2-parameter family of symmetric integrals of 
degree 16. 

For example, take d 3 — d 4 — 1 to obtain 

z 4 

(7.7) Rx(z) = ( z ie + z i4 + z u + 115z io + 20z 8 + 115z 6 + z 4 + z 2 + l) 2 ' 
The main theorem yields 

1024z 4 + 2304z 6 + 1792z 8 + 560z 10 + 60z 12 + z 14 



(7.8) R 2 {z) - 



2 7 (16z 8 + 36z 6 + 27z 4 + 36z 2 + 16) 5 
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SO that 

/•oo /*oo 

/ Ri(z) dz = / R,2(z)dz. 
Jo Jo 

Letting f[n] := N n ,g[l, n, 27/32, 9/4] we obtain 

R x {z)dz = 2- 15 (/[0]+60/[l] + 1584/[2] + 4096/[3]) 



o 

and conclude that 

/•OO 







(z 16 + z 14 + z 12 + 115z 10 + 20z 8 + 115z 6 + z 4 + z 2 

(149288517+ 12947003v / T31> 
1124663296V54925 + 4798V131 ' 



Example 12. We classify the symmetric denominators of degree 32 that yield 
computable integrals. These functions depend on 8 parameters 



(7.9) 



D 8 (d u --- ,d 8 ;z) = ]Td 9 ^(z 2fc + z 2 ( 16 - fe )) 



k=0 



and the main theorem expresses the integral in terms of Eg. The conditions for Eg 
to be symmetric yield 

d l = -3441 + 35d 5 + 64d 6 - 312d 7 - 3264d 8 

d 2 = 34720 - 56d 5 - 110d 6 + 560c2 7 + 4565d 8 

d 3 = -3472 + 28d 5 + 64d 6 - 329d 7 - 2240d 8 

d 4 = 4960 - 8d 5 - 19rf 6 + 8Qd? + 938<i 8 

and the symmetric Eg is 

E s (d 5 ,d 6 ,d 7 , d a ;z) = 32768(1 + z 16 ) + (131072 + 8192d 8 )(z 2 + z 14 )+ 

(212992 + 2048d 7 + 28672d 8 )(z 4 + z 12 ) + (180224 + 512d 6 + 6144d 7 + 39424d 8 )(z 6 + z 10 )- 

+ (84480 + 128rf 5 + 1280d 6 + 6912d 7 + 26880d 8 )z 8 . 
The symmetry of Eg now determines d^ 1 d§ in terms of c? 7 , dg and we obtain 



(7.10) 



The function (z 2 + l) 16 is a symmetric polynomial of degree 32 and yields a par- 
ticular solution to (|7.10p . 
As before let 



'di 






" 63475 " 




9166 " 




' 54640 " 


d 2 






-100800 




-14392 




-86645 


d 3 




-31 


47936 




6895 




41664 


di 




-13664 


+ 


-1964 


d 7 + 


-11471 


d 5 






2220 




322 




2000 


_d 6 _ 






-224 




-28 




-189 



= d 3 - 



16 



for 2 < j < 8 and e\ = d\ — \ 
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Then 



(7.11) 



ei 




9166 " 




" 54640 " 


e2 




-14392 




-86645 


e.3 




6895 




41664 


e 4 




-1964 


e 7 + 


-11471 


65 




322 




2000 


_ e 6. 




-28 




-189 



and as in the case of degree 16 we can compute a 2-parameter family of symmetric 
integrals of degree 32. 



Appendix A. Two binomial sums 

The closed-form evaluation of sums involving binomials coefficients can be ob- 
tained by traditional analytical techniques or by using the powerful WZ-method 
as described in [16]. We discuss two sums used to simplify expressions in later 
sections, presenting one proof in each style. 

Lemma A.l. Let k, N be positive integers with k < N. Then 

gcvx;o - rry-* 

Proof. Multiply the left hand side of (|A.1[) by x k and sum over k to produce 



k=Q j=0 J / \ / j =Q \ J / k=Q 



x k 



3=0 

(l + ^x~TT) 2N+1 - (l - V^TT) 2JV+1 



2^/x~TT 



,N-k x k 



k=0 v 

In order to justify the last step we start with the well known result 



A- 



ttM^f) ■ ?cv> 

(see WILF [19], page 54). Letting x = -Ay and i = 2N + 1 in (|A.2j) gives 
(i-V^+T) 2N+1 fv_ k+1 [2N + l + 2k\ {N+1+k) 2N+1+k . 

similarly x — —Ay and i = — 2N — 1 yields 



(l + V^ + T) 2 ^ 1 = g ( _ i)fe+1 /2A+l + 2fc 

v A:=2JV+1 v 



-(.N+l+k) 2N+1+H 
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Thus 

2iV 



(l + v^ + Tr+i-g-vCT) 2 ^ 1 ^ f-(2N + l-2k)\ 

2yGT+T ~ ^ M fc j 4 



fc=0 



Lemma A.2. Let N e N. Then 



fc=0 



□ 



j=o V J / j=0 v J 7 



Proof. The coefficient of z on the left hand side is 

N-k 

E 



2N+ l\ (N — j 
2j ){ k 



and the corresponding coefficient on the right hand side is ( 2N , k )A N k . The result 



then follows from Lemma [A. II □ 
N e N with k < N. Then 
/2N\/N-j\ 2 2k ~ l N fk + N-V 



Lemma A.3. Let fc, N € N with k < N. Then 

k 

E \ 2/ M :v - /,■/ /.• \ -V • • /. 

Proof. This lemma could be proven in the same style as Lemma IA.ll Instead we 
use the WZ- method as explained in [16]. Indeed, let 

F{k-j) - 



N2^( k + N -r 
and define, with the package EKHAD, the function 

m - 1) 



G(k;j) = F(k;j)x 



2(N + k)(k-j + l) 

Then F(k;j) — F(k + 1; j) = G(k;j + 1) — G(k;j), and summing over j we see that 
the sum of F{k\ j) over j is independent of fc. The case fc = N produces 1 as the 
common value. □ 

Lemma A. 4. Let p € N, d\, cfe, • • • } dp be parameters, and define dp+i := 1. Then 

(A.3) t^iV* Ef 2P 2 "/ fc ) (1 + z2) ^ = 



fc=0 j=o 
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Proof. For fixed < i < p — 1 the coefficient of z 2t on the right hand side of (|A.3j) 
is 

^ r=p-i+l \ -P 1 / 

and for i = p we have [RHS] (2p) = 1 + Y^j=i dj- Similarly, for the left hand side 
of (TO). 

P+i 

[LHS](2i) = d 

r=p-\-l—i 

It is easy to check that the coefficients of z 2p match. It suffices to show that for 
each i such that < i < p — 1 and for each r such that p+l-j<r<p+lwe 
have 

^/2r-2\/ r-l-j \ 2 2 (P" l )~ 1 fr+p-i-2 

p-" Q \ 2 j / \r - 1 - p + i) ~ p - i ' \r - p + i - 1, 

This follows from Lemma IA.1I with k = p — i and iV = r — 1. □ 

The suggestions of the referrees and the editor are gratefully acknowledged. 

References 




[12: 

[13; 
[14; 

[15 



BOROS, G. - MOLL, V.: A criterion for unimodality. Elec. Jour, of Combinatorics, 6 (1999), 
#R10. 

BOROS, G. - MOLL, V.: An integral hidden in Gradshteyn and Ryzhik. Jour. Comp. Appl. 
Math. 106, 361-368, 1999. 

BOROS, G. - MOLL, V.: A rational Landen transformation. Contemporary Mathematics 
251, 83-91, 2000. 

BORWEIN, J. - BORWEIN, P.: Pi and the AGM. Canadian Mathematical Society. Wiley- 
Interscience Publication. 

BRONSTEIN, M.: Symbolic Integration I. Transcendental functions. Algorithms and Com- 
putation in Mathematics, 1. Springer- Verlag, 1997. 

GEDDES, K. - CZAPOR, S.R. - LABAHN, G.: Algorithms for Computer Algebra. Kluwer, 
Dordrecht. The Netherlands, 1992. 

GAUSS, K.F.: Arithmetische Geometrisches Mittel, 1799. In Werke, 3, 361-432. Konigliche 
Gesellschaft der Wissenschaft, Gottingen. Reprinted by Olms, Hildcscheim, 1981. 
GRADSHTEYN, I.S. - RYZHIK, I.M.: Table 0} Integrals, Series and Products. Fifth Edition, 
ed. Alan Jeffrey. Academic Press, 1994. 

HARDY, G.H.: The Integration of Functions of a Single Variable. Cambridge Tracts in 
Mathematics and Mathematical Physics, 2, Second Edition, Cambridge University Press, 
1958. 

HERMITE, C: Sur I'integration des fractions rationelles. Nouvelles Annales de Mathema- 
tiqucs ( 2 Gmc serie) 11, 145-148, 1872. 

HOROWITZ, E.: Algorithms for partial fraction decomposition and rational function inte- 
gration. Proc. of SYMSAM'71, ACM Press, 441-457, 1971. 

LAZARD, D. - RIOBOO, R,.: Integration of Rational Functions: Rational Computation of 
the Logarithmic Part. Journal of Symbolic Computation 9, 113-116, 1990. 
MCKEAN, H. - MOLL, V.: Elliptic Curves: Function Theory, Geometry, Arithmetic. Cam- 
bridge University Press, 1997. 

NEWMAN, D.: A simplified version of the fast algorithm of Brent and Salamin. Math. 
Comp. 44, 207-210, 1985. 

OSTROGRADSKY, M.W.: De I'integration des fractions rationelles. Bulletin de la Classe 
Physico-Mathematiques de l'Acadcmic Imperieriale des Sciences de St. Petersbourgh, IV, 
145-167, 286-300. 1845. 



20 



GEORGE BOROS AND VICTOR H. MOLL 



[16] PETKOVSEK, M. - WILF, H.S. - ZEILBERGER, D.: A=B. A. K. Peters, Wellesley, Mas- 
sachusetts. 1996. 

[17] ROTHSTEIN, M.: A new algorithm for the integration of Exponential and Logarithmic 
Functions, Proc. of the 1977 MACSYMA Users Conference, NASA Pub., CP-2012, 263-274. 

[18] TRACER, B.M.: Algebraic factoring and rational function integration. Proc. SYMSAC 76, 
219-226. 

[19] WILF, H.S.: generatingfunctionology. Academic Press, 1990. 

Department of Mathematics, Xavier. University, New Orleans, Louisiana 70125 
E-mail address: gborosSxula.math.edu 

Department of Mathematics, Tulane University, New Orleans, LA 70118 
E-mail address: vhmamath . tulane . edu 



